ratio.se <- function(num, den, num.var, den.var, nd.cov) {
  vv <- (1/den^2) * num.var + ((num^2)/(den^4)) * den.var - 2 * (num/(den^3)) * nd.cov
  return(sqrt(vv))
}

joint.iv <- function(y, d1, d2, z1, z2) {
  i11 <- z1 == 1 & z2 == 1
  i10 <- z1 == 1 & z2 == 0
  i01 <- z1 == 0 & z2 == 1
  i00 <- z1 == 0 & z2 == 0
  n11 <- sum(i11)
  n10 <- sum(i10)
  n01 <- sum(i01)
  n00 <- sum(i00)
  N <- length(y)
  d11 <- d1 * d2
  d10 <- d1 * (1 - d2)
  d01 <- (1 - d1) * d2
  d00 <- (1 - d1) * (1 - d2)

  rho <- rep(NA, times = 9)
  rho.var <- rep(NA, times = 9)
  names(rho) <- c("cc", "ca", "cn", "ac", "aa", "an", "nc", "na", "nn")
  names(rho.var) <- c("cc", "ca", "cn", "ac", "aa", "an", "nc", "na", "nn")
  rho["cc"] <- mean(d11[i11]) - mean(d11[i10]) - mean(d11[i01]) + mean(d11[i00])
  rho["ca"] <- mean(d11[i10]) - mean(d11[i00])
  rho["cn"] <- mean(d00[i01]) - mean(d00[i11])
  rho["ac"] <- mean(d11[i01]) - mean(d11[i00])
  rho["aa"] <- mean(d11[i00])
  rho["an"] <- mean(d10[i01])
  rho["nc"] <- mean(d00[i10]) - mean(d00[i11])
  rho["na"] <- mean(d01[i10])
  rho["nn"] <- mean(d00[i11])
  rho.var["cc"] <- (N/(N-1)) * (var(d11[i11])/n11 + var(d11[i10])/n10 + var(d11[i01])/n01 + var(d11[i00])/n00)
  rho.var["ca"] <- (N/(N-1)) * (var(d11[i10])/n10 + var(d11[i00])/n00)
  rho.var["cn"] <- (N/(N-1)) * (var(d00[i01])/n01 + var(d00[i11])/n11)
  rho.var["ac"] <- (N/(N-1)) * (var(d11[i01])/n01 + var(d11[i00])/n00)
  rho.var["aa"] <- (N/(N-1)) * (var(d11[i00])/n00)
  rho.var["an"] <- (N/(N-1)) * (var(d10[i01])/n01)
  rho.var["nc"] <- (N/(N-1)) * (var(d00[i10])/n10 + var(d00[i11])/n11)
  rho.var["na"] <- (N/(N-1)) * (var(d01[i10])/n10)
  rho.var["nn"] <- (N/(N-1)) * (var(d00[i11])/n11)

  tau <- rep(NA, times = 12)
  names(tau) <- c("cc.10", "cc.11", "cc.20", "cc.21", "cn", "ca", "nc", "ac", "laje", "laie", "cc.cn", "cc.nc")
  tau.se <- tau

  ## cc.10
  yp <- y * (1-d2)
  num <- mean(yp[i10]) - mean(yp[i11]) - mean(yp[i00]) + mean(yp[i01])
  num.var <- var(yp[i10])/n10 + var(yp[i11])/n11 + var(yp[i00])/n00 + var(yp[i01])/n01
  num.var <- (N/(N - 1)) * num.var
  den <- rho["cc"]
  den.var <- rho.var["cc"]
  nd.cov <- -cov(yp[i10], d11[i10])/n10 - cov(yp[i11],d11[i11])/n11 -
    cov(yp[i00], d11[i00])/n00 - cov(yp[i01], d11[i01])/n01
  nd.cov <- (N/(N - 1)) * nd.cov
  tau["cc.10"] <- num/den
  tau.se["cc.10"] <- ratio.se(num, den, num.var, den.var, nd.cov)

  ## cc.11
  yp <- y * d2
  num <- mean(yp[i11]) - mean(yp[i10]) - mean(yp[i01]) + mean(yp[i00])
  num.var <- var(yp[i11])/n11 + var(yp[i10])/n10 + var(yp[i01])/n01 + var(yp[i00])/n00
  num.var <- (N/(N - 1)) * num.var
  den <- rho["cc"]
  den.var <- rho.var["cc"]
  nd.cov <- cov(yp[i11], d11[i11])/n11 + cov(yp[i10],d11[i10])/n10 +
    cov(yp[i01], d11[i01])/n01 + cov(yp[i00], d11[i00])/n00
  nd.cov <- (N/(N - 1)) * nd.cov
  tau["cc.11"] <- num/den
  tau.se["cc.11"] <- ratio.se(num, den, num.var, den.var, nd.cov)

  ## cc.20
  yp <- y * (1-d1)
  num <- mean(yp[i01]) - mean(yp[i11]) - mean(yp[i00]) + mean(yp[i10])
  num.var <- var(yp[i01])/n01 + var(yp[i11])/n11 + var(yp[i00])/n00 + var(yp[i10])/n10
  num.var <- (N/(N - 1)) * num.var
  den <- rho["cc"]
  den.var <- rho.var["cc"]
  nd.cov <- -cov(yp[i10], d11[i10])/n10 - cov(yp[i11],d11[i11])/n11 -
    cov(yp[i00], d11[i00])/n00 - cov(yp[i01], d11[i01])/n01
  nd.cov <- (N/(N - 1)) * nd.cov
  tau["cc.20"] <- num/den
  tau.se["cc.20"] <- ratio.se(num, den, num.var, den.var, nd.cov)

  ## cc.21
  yp <- y * d1
  num <- mean(yp[i11]) - mean(yp[i10]) - mean(yp[i01]) + mean(yp[i00])
  num.var <- var(yp[i11])/n11 + var(yp[i10])/n10 + var(yp[i01])/n01 + var(yp[i00])/n00
  num.var <- (N/(N - 1)) * num.var
  den <- rho["cc"]
  den.var <- rho.var["cc"]
  nd.cov <- cov(yp[i11], d11[i11])/n11 + cov(yp[i10],d11[i10])/n10 +
    cov(yp[i01], d11[i01])/n01 + cov(yp[i00], d11[i00])/n00
  nd.cov <- (N/(N - 1)) * nd.cov
  tau["cc.21"] <- num/den
  tau.se["cc.21"] <- ratio.se(num, den, num.var, den.var, nd.cov)

  ## cn
  yp <- y * (1 - d2)
  num <- mean(yp[i11]) - mean(yp[i01])
  num.var <- var(yp[i11])/n11 + var(yp[i01])/n01
  num.var <- (N/(N - 1)) * num.var
  den <- mean(d00[i01]) - mean(d00[i11])
  den.var <- (N/(N-1)) * (var(d00[i01])/n01 + var(d00[i11])/n11)
  nd.cov <- -cov(yp[i11], d00[i11])/n11 - cov(yp[i01],d00[i01])/n01
  nd.cov <- (N/(N - 1)) * nd.cov
  tau["cn"] <- num/den
  tau.se["cn"] <- ratio.se(num, den, num.var, den.var, nd.cov)

  ## ca
  yp <- y * d2
  num <- mean(yp[i10]) - mean(yp[i00])
  num.var <- var(yp[i10])/n10 + var(yp[i00])/n00
  num.var <- (N/(N - 1)) * num.var
  den <- mean(d11[i10]) - mean(d11[i00])
  den.var <- (N/(N-1)) * (var(d11[i10])/n10 + var(d11[i00])/n00)
  nd.cov <- cov(yp[i10], d11[i10])/n10 + cov(yp[i00],d11[i00])/n00
  nd.cov <- (N/(N - 1)) * nd.cov
  tau["ca"] <- num/den
  tau.se["ca"] <- ratio.se(num, den, num.var, den.var, nd.cov)

  ## nc
  yp <- y * (1 - d1)
  num <- mean(yp[i11]) - mean(yp[i10])
  num.var <- var(yp[i11])/n11 + var(yp[i10])/n10
  num.var <- (N/(N - 1)) * num.var
  den <- mean(d01[i11]) - mean(d01[i10])
  den.var <- (N/(N-1)) * (var(d01[i11])/n11 + var(d11[i10])/n10)
  nd.cov <- cov(yp[i11], d11[i11])/n11 + cov(yp[i10],d11[i10])/n10
  nd.cov <- (N/(N - 1)) * nd.cov
  tau["nc"] <- num/den
  tau.se["nc"] <- ratio.se(num, den, num.var, den.var, nd.cov)

  ## ac
  yp <- y * d1
  num <- mean(yp[i01]) - mean(yp[i00])
  num.var <- var(yp[i01])/n01 + var(yp[i00])/n00
  num.var <- (N/(N - 1)) * num.var
  den <- mean(d11[i01]) - mean(d11[i00])
  den.var <- (N/(N-1)) * (var(d11[i01])/n01 + var(d11[i00])/n00)
  nd.cov <- cov(yp[i01], d11[i01])/n01 + cov(yp[i00],d11[i00])/n00
  nd.cov <- (N/(N - 1)) * nd.cov
  tau["ac"] <- num/den
  tau.se["ac"] <- ratio.se(num, den, num.var, den.var, nd.cov)

  ## laje
  yp <- y * (d11 - d00)
  num <- mean(yp[i11]) - mean(yp[i10]) - mean(yp[i01])  + mean(yp[i00])
  num.var <- var(yp[i11])/n11 + var(yp[i10])/n10 + var(yp[i01])/n01 + var(yp[i00])/n00
  num.var <- (N/(N - 1)) * num.var
  den <- rho["cc"]
  den.var <- rho.var["cc"]
  nd.cov <- cov(yp[i11], d11[i11])/n11 + cov(yp[i10], d11[i10])/n10 + cov(yp[i01], d11[i01])/n01 + cov(yp[i00],d11[i00])/n00
  nd.cov <- (N/(N - 1)) * nd.cov
  tau["laje"] <- num/den
  tau.se["laje"] <- ratio.se(num, den, num.var, den.var, nd.cov)

  ## laie
  yp <- y
  num <- mean(yp[i11]) - mean(yp[i10]) - mean(yp[i01])  + mean(yp[i00])
  num.var <- var(yp[i11])/n11 + var(yp[i10])/n10 + var(yp[i01])/n01 + var(yp[i00])/n00
  num.var <- (N/(N - 1)) * num.var
  den <- rho["cc"]
  den.var <- rho.var["cc"]
  nd.cov <- cov(yp[i11], d11[i11])/n11 + cov(yp[i10], d11[i10])/n10 + cov(yp[i01], d11[i01])/n01 + cov(yp[i00],d11[i00])/n00
  nd.cov <- (N/(N - 1)) * nd.cov
  tau["laie"] <- num/den
  tau.se["laie"] <- ratio.se(num, den, num.var, den.var, nd.cov)

  ## cc.cn
  yp <- y * (1 - d2)
  num <- mean(yp[i10]) - mean(yp[i00])
  num.var <- var(yp[i10])/n10 + var(yp[i00])/n00
  num.var <- (N/(N - 1)) * num.var
  den <- mean(d10[i10]) - mean(d10[i00])
  den.var <- (N/(N-1)) * (var(d10[i10])/n10 + var(d10[i00])/n00)
  nd.cov <- cov(yp[i10], d10[i10])/n10 + cov(yp[i00],d10[i00])/n00
  nd.cov <- (N/(N - 1)) * nd.cov
  tau["cc.cn"] <- num/den
  tau.se["cc.cn"] <- ratio.se(num, den, num.var, den.var, nd.cov)

  ## cc.nc
  yp <- y * (1 - d1)
  num <- mean(yp[i01]) - mean(yp[i00])
  num.var <- var(yp[i01])/n01 + var(yp[i00])/n00
  num.var <- (N/(N - 1)) * num.var
  den <- mean(d01[i01]) - mean(d10[i00])
  den.var <- (N/(N-1)) * (var(d01[i01])/n01 + var(d10[i00])/n00)
  nd.cov <- cov(yp[i01], d10[i01])/n01 + cov(yp[i00],d10[i00])/n00
  nd.cov <- (N/(N - 1)) * nd.cov
  tau["cc.nc"] <- num/den
  tau.se["cc.nc"] <- ratio.se(num, den, num.var, den.var, nd.cov)


  return(list(tau=tau, se = tau.se, rho = rho, rho.se = sqrt(rho.var)))
}
